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Abstract. We prove a necessary and sufficient criterion for the exponential 
I i stability of periodic solutions of delay differential equations with large delay. 

We show that for sufficiently large delay the Floquet spectrum near criticality is 
^*\ characterized by a set of curves, which we call asymptotic continuous spectrum, 

that is independent on the delay. 

•5 

1. Introduction. Delay-differential equations (DDEs) are similar to ordinary dif- 
ferential equations (ODEs) except that the right-hand side may depend on the past, 
i— i For example, they could be of the form 

x(t) = f(x(t),x(t-r)) (1) 

where x(t) is a vector in K™ and the delay r > decides how far one looks into the 
past. When studying DDEs as dynamical systems one notices that equilibria do not 
y—i depend on the delay r. More precisely, their location and number is independent 

t-H of t. However, their stability changes significantly when one varies r, an effect 

that is well known and of practical importance in engineering and control [6, 12]. 
Exponential stability is given by the spectrum of the linearization of the DDE in 
its equilibrium. This spectrum, in turn, can be expressed as roots of an analytic 
function (a polynomial of exponentials A t— > exp(— At)). Lichtner et al [5] classified 
^ rigorously which types of limits this spectrum can have as r approaches infinity. 

I^J Roughly speaking, for sufficiently large r all except maximally n eigenvalues form 

bands near the imaginary axis. After rescaling their real part by 1/r one finds 
that these bands converge to curves, called asymptotic continuous spectrum. They 
are given as root curves of parametric polynomials, and are, thus, much easier to 
compute than the eigenvalues of the singularly perturbed large-delay problem. Of 
practical relevance are then stability criteria based entirely on the asymptotic spectra 
that guarantee the stability of an equilibrium for sufficiently large delays r. 

This paper gives a similar result for periodic orbits of (1). In contrast to equi- 
libria, periodic orbits change as the delay r varies, such that the statement about 
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Figure 1. Illustration of stability for the example (2)-(3) (com- 
puted with DDE-Biftool [1, 7, 8]): (a) bifurcation diagram of peri- 
odic orbits in the (r, T)-plane; (b) phase portraits of periodic orbits 
along branch; (c,d) images of the branch under transformation 
t H> NT(t) + r, illustrating overlaps and coexistence; (e) domi- 
nant Floquet exponents of selected periodic orbit (grey square in 
panel (a)); (f) dominant Floquet exponents after rescaling of real 
parts and curves of asymptotic continuous spectrum. Parameters: 
a = -0.10779 for (a-d), r = 1.081, N = 25 for (c,d). 



"independence" of the delay (or, rather, re-appearance) has to be formulated more 
carefully. Let us look at a two-dimensional example to illustrate the observation 
made by Yanchuk & Pcrlikowski [16]: 

x(t) = ax(t) - 2iry(t) - x{t) [x{tf + y{tf] (2) 
y{t) = 2nx{t) + ay(t) - y(t) [x{tf + y{t) 2 ] + y(t - r) (3) 

where we fix a ~ —0.1 and vary r. System (2)-(3) consists of the normal form for 
the supercritical Hopf bifurcation with an additional delayed term y(t — r) in the 
second equation, which breaks the rotational symmetry of the instantaneous terms. 
Numerically, one can observe that the system has a family of periodic orbits for delays 
between r w 0.7 and r w 1.3. The period and the phase portrait projected onto the 
(x,y)-plane of these orbits are shown in Figure 1(a) and (b). The family repeats, the 
orbits keeping their shape, for every integer N if we change the horizontal axis in 
Figure 1(a) to NT(t) + r. The transformation NT(t) + t is not just a parallel shift, 
since the dependence T(r) is, generically, nontrivial (in the example dT/dr > 0). 
This leads to an overlapping of the families for large r as shown in Figure 1(c) and 
(d). 
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Specifically, let us consider a regular periodic orbit Xq, existing for a fixed value To 
of r (regularity means that the unit Floquet multiplier of xq is simple, see Section 2) . 
Then xq persists for r in a small neighborhood of To such that we have a branch of 
periodic orbits depending locally on the parameter r close to To. Let U be a small 
neighborhood of this orbit along the branch. Along this branch the period T of the 
orbit is a smooth function T{t) of t. In Figure 1(a) an example orbit Xo is indicated 
by a grey square at r = 1.081 and its neighborhood U is highlighted by a slight 
thickening of the line. If we assume that T'(tq) ^ (which is a genericity condition) 
then T(t) will have a slope uniformly bounded away from zero for all r in the 
neighborhood corresponding to U. This implies that the images of the neighborhood 
U are stretched proportionally to N under the transformation r i-> NT{t) + t. Thus, 
for any given large r a large number of periodic orbits from U coexist, and the 
number of coexisting orbits is proportional to r (see the overlapping images of IA in 
Figure 1(c) and (d)). 

The next question to ask is: which of those many coexisting periodic orbits for 
large delays are dynamically stable? Is it possible to derive sharp stability criteria for 
the large-delay orbits based on quantities independent of the delay? More precisely, 
what is the stability of a given periodic orbit xq (such as the one indicated by a grey 
square in Figure 1(a) for delays t = NT(t ) + t as N goes to infinity (this would 
be the sequence of orbits indicated by grey squares in Figure l(c,d)). 

Our example suggests that the stability for t = NT(tq) + To where N — > oo can 
indeed be determined from spectral properties of the periodic orbit at the small 
delay To- Figure 1(e) shows the spectrum of the periodic orbit xq highlighted by a 
grey square in Figure 1(a) for N — 25. One can see that, first, its Floquet exponents 
form bands, and, second, there is a large number of Floquet exponents very close to 
the imaginary axis (note the scale of the horizontal axis in Figure 1(e)). Figure 1(f) 
illustrates one of the results of this paper: after rescaling the real part, Floquet 
exponents converge to curves for N — > oo. These curves, the asymptotic continuous 
spectrum are computable by solving regular periodic boundary value problems 
parametrized by w, the vertical axis in Figure 1(f). Since the original nonlinear DDE 
is autonomous, one of the curves of the asymptotic continuous spectrum touches the 
imaginary axis. Similar to the equilibrium case, we establish that the asymptotic 
continuous spectrum (together with the strongly unstable spectrum, see Section 2) 
determines the stability of the periodic orbit xq for sufficiently large N. 

We have colored the branch in Figure 1(a) already according to the conclusions 
from the asymptotic spectra. The part that is displayed as a black curve in Figure 1(a) 
is the part of the branch that will be exponentially stable as N tends to infinity, 
whereas the grey part will be exponentially unstable, having a large number of 
weakly unstable Floquet exponents. 

In short, the observation by Yanchuk & Perlikowski [1G] implies that, if we find 
a periodic orbit xo of the DDE (1) for a fixed small delay To, then (under some 
genericity conditions) the DDE (1) has a large number of similar orbits coexisiting 
for any sufficiently large delay r. This paper establishes a sharp critierion (again 
under some genericity conditions) determining if all of these coexisting periodic 
orbits are dynamically stable. It does so by providing formulas for the Floquet 
exponents of the periodic orbit with delay r = NT(t q ) + t that are independent of 
N but valid asymptotically for large N. For the example DDE (2)-(3), our results 
imply that the system has a large number of coexisting stable periodic orbits for 
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every sufficiently large delay r (as suggested by the Figures l(c,d) where the black 
parts of the branches are stable). 

Section 2 gives a non-technical overview of the results gradually developed and 
proven in the later sections. One central part of our paper is the construction of a 
characteristic function 

H h> h (jm, cx P (-(NT{t ) + tq)h)) , 

the roots of which are the Floquet exponents of the periodic orbit for delay To + 
JVT(to), and for which we can study the limit N — > oo. This construction is given 
in Section 3. The existence of this function h permits us to follow the approach 
from [5] and to extend their techniques to the case of periodic solutions. In the 
following sections 4 and 5 we describe two parts of the Floquet spectrum that 
show a different scaling behavior for large N (and, thus, r). The strongly unstable 
spectrum, converging to a finite number of asymptotic Floquet exponents which 
are determined by the instantaneous terms, is investigated in Section 4. Then, in 
Section 5, we analyze the Floquet exponents given by the asymptotic continuous 
spectrum shown in Fig 1(f). Based on these results, we can then prove in Section 6 
our main result, a criterion about asymptotic stability based on the location of the 
asymptotic continuous and strongly unstable spectrum. 

In contrast to the case of spectra at equilibria, where the asymptotic continuous 
and strongly unstable spectrum in many cases can be calculated explicitly (see [5]), 
the corresponding parts of the Floquet spectrum of a periodic orbit xq, can typically 
only be computed numerically. This limitation is not specific to our results but 
applies equally to most stability results based on the Floquet spectrum for ODEs. 
Even the periodic orbit itself is in most cases only computable with numerical 
methods. An exception is the case, where the periodic orbit is at the same time a 
symmetry orbit of the system. In this case, examples of asymptotic continuous and 
strongly unstable Floquet spectrum have been calculated explicitly in [15, 17]. 

Conversely, the results presented in this paper provide an approach to approximate 
Floquet spectra numerically for large delays r. If one uses numerical methods on 
problems with large delays, one faces the difficulty that the size of the matrix arising 
in the discretized eigenvalue problem grows not only with the desired accuracy 
(which is natural) but also with r, even if the period of the orbit remains bounded. 
This is the case for the numerical methods used in DDE-Biftool [7]. This increase 
is to be expected because the number of Floquet exponents close to the imaginary 
axis increases with r (see Theorem 6). In contrast to this, the asymptotic spectra 
can be computed with the same numerical method as in DDE-Biftool and a matrix 
size that is uniformly bounded for large r. This paper does not discuss the details 
of the numerical computation of asymptotic spectra. However, our construction of 
the characteristic function h is uniform with respect to r and could in principle be 
implemented numerically. In practice, it is better to apply the same analysis as is 
done in this paper to the large discretized eigenvalue problem. 

In Fig. 1(f) we demonstrated the good agreement between original Floquet 
exponents and asymptotic spectra already in the situation of a moderately large 
delay, where the original Floquet exponents are still computable. Moreover, our 
stability criterion, which can be reliably verfied in this way, is valid independently 
on the actual value of the large delay. 

2. Basic concepts and overview of the results. 
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2.1. Periodic orbits, stability, and Floquet exponents. Let x* be a periodic 
orbit of the n-dimensional autonomous nonlinear delay differential equation (DDE) 

x(t)=f(x(t),x(t- T )), (4) 

that is, x*(t) satisfies (4) for all times t and has period T: x*(t) — x*(t + T) for all 
t € R. Without loss of generality we may assume that T = 1 (this can be achieved 
by a rescaling of time and the delay r). If the period T is equal to 1 then x* is also 
a periodic orbit of 

x(t) = f(x(t),x(t-r-N)), (5) 

where N is a natural number and r € [0, 1). We denote the restriction of the periodic 
function x* : M. n- R n to the interval [— t — TV, 0] also by x* such that x* is an 
element of C([-r - AT, 0];R n ). 

We are concerned with the question whether the periodic orbit x* is stable or 
unstable for sufficiently large N in the following sense: 

Definition 1 (Exponential orbital stability and instability). 

Let X(t; •) be the semiflow on C([-r - N, 0]; R n ) induced by DDE (5). The periodic 
orbit is called exponentially orbitally stable if there exists a decay rate 7 > 
such that all initial history segments Xq in a neighborhood of x* satisfy 

\\X(t; x a ) - X(t + t Q ; x*)\\oo < C exp(-7i)||x - x*||oo 

for some time shift to and some constant C > 0. 

Similarly, x* is called exponentially unstable if there exists a growth rate 
7 > 0, a neighborhood A/" of x* and a constant C > such that one can find initial 
history segments xo ^ x* arbitrarily close to x* that satisfy 

\\X(n:x Q ) - XfcxJWoo > Cexp(7n)||x - > 

for all n e N as long as X(n; xq) stays in the neighborhood M . 

This is the standard definition for stability of periodic orbits used also for ODEs 
except that the phase space is C([— r — N, 0];M. n ) instead of W 1 . The notation || • 
refers to the usual maximum norm in C([— r — N, 0];M n ). 

Textbook theory of delay equations reduces the stability problem to the problem of 
finding eigenvalues of the linear map M N : C([-r - N, 0]; C") -> C([-r - N, 0]; C"), 
which is given as the time-1 map of the linear DDE 

±(t) = A(t)x(t) + B(t)x(t-T-N), (6) 

where the time-dependent n x n-matrices A(t) € M™ xn and S(t) G R nx ™ are the par- 
tial derivatives of / in x*: A(<) = dif(x*(t), x»(t — r)) and f?(t) = d2f(x*(t),x*(t — 
r)) [2]. If the right-hand side / of the nonlinear problem (4) is smooth in its 
arguments then the matrices A and B are also smooth periodic function of time t in 
the interval [0,1]. 

Since the time-1 map (Mjy) N+1 is compact the spectral theory for compact 
operators and the polynomial spectral mapping theorem imply that the spectrum 
cr(-Mjv) consists of a sequence of eigenvalues of finite multiplicity accumulating only 
at zero (zero is the only element of ct(Mjv) that is not in the point spectrum). Also, 
A = 1 is always an eigenvalue of Mm because x*(i) satisfies (6) and has period 1. 

The periodic orbit x* is exponentially {orbitally) stable if and only if 

(1) the eigenvalue 1 of Mm is algebraically simple, and 

(2) all other eigenvalues of Mm have modulus less than 1. 
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Similarly, it is exponentially unstable if at least one eigenvalue has modulus greater 
than 1. 

Thus, the stability of iterations of Mm is determined by its eigenvalues. We also 
use the term orbitally stable for the map Mjv, meaning that Mm satisfies both of 
the above conditions. 

The state space of DDE (6) is the function space C([— r — N, 0];R"). Thus, initial 
value problems for (6) require specifying an infinite-dimensional initial condition. 
Similarly, one expects that a boundary value problem for a DDE requires the 
specification of infinitely many boundary conditions. However, periodic boundary 
value problems are easier to formulate: for example, a solution of the periodic 
boundary value problem for the general DDE (4) for period T = 1 is simply a 
function x £ C l ([—1, 0];K n ) satisfying 



where the notation (t — T) mod [_ 1 ] stands for t — r + fcifi — t £ [—k — 1, —k) 
and k £ Z is the integer part of r — t. Since x is continuous and satisfies the 
periodicity condition (8) it can be extended continuously to a continuous function 
on the whole real line by defining x(s) = x{s mo d[-i,o])- Consequently, the right-hand 
side of (7) is continuous for alH £ [—1,0], which guarantees that x can really satisfy 
the differential equation pointwise and is an element of C . The solution x then 
automatically satisfies x(— 1) = x(Q), and, thus, by induction is as smooth as the 
right-hand side /. In this respect, periodic boundary value problems for DDEs are 
similar to boundary value problems for ODEs. In Section 3 we will reduce linear 
periodic boundary value problems for DDEs to low-dimensional linear systems of 
algebraic equations. 

Definition 2. We call /i a Floquet exponent of Mm, and write [i £ £jv, if exp(/i) is 
an eigenvalue of Mm- 

Floquet exponents of Mm can be found as those complex numbers fi for which 
the periodic boundary value problem 



has a nontrivial solution y £ C 1 ([— 1, 0], C" ) [2]. Note that in (9) we use the delay 
r £ [0,1) to calculate the Floquet exponents of Mm, i.e. for a periodic orbit of 
system (6) with delay t + N. Only the factor exp(— Nfj,) in front of B(t) accounts 
for the large delay whereas we have just extended y periodically for arguments less 
than —1. If (9) has a non-trivial solution y(t) for /i £ C then it has the non-trivial 
solution y(t) exp(2irikt) for /j, + 2kir\ for any integer k. Hence, we choose the Floquet 
exponent [i such that its imaginary part is between [— 7r, 7r). 

2.2. Asymptotic spectra for TV — > oo. The set of Floquet exponents, Eat, forms 
a discrete subset of the complex plane, the point spectrum of exponents, which 
depends on N. In order to describe in which form Y>m has a limit for N — > oo we 
introduce two asymptotic spectra, which are also subsets of the complex plane. The 
notation follows [16]. 



= f(x{t), x(t - r) mod [-i,o]) 
x(-l) = x(0) 



(7) 
(8) 



ij{t) - [A(t) - fxl]y{t) + exp(-(iV + r) M ) B(t) y((t - r) mod[ _ 1>0] ) 
1/(0) = !/(-!) 



(9) 
(10) 
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Definition 3 (Instantaneous and strongly unstable spectrum). The set of all 
/i G C for which the linear ODE boundary value problem on [—1, 0] 



has a non-trivial solution y E C 1 ([— 1, 0]; C") is called the instantaneous spectrum. 
The subset A+ C Yja of those fi with positive real part is called the strongly 
unstable asymptotic spectrum. 

The instantaneous spectrum Tia contains exactly n elements with imaginary part 
in [— 7T, 7t), counting algebraic multiplicity. We note that Tia and A+ do not depend 
on N but only on A. One result of our paper is that all Floquet exponents of Mjy 
with a real part that is positive uniformly in N converge to elements of the strongly 
unstable spectrum A+. 

Yanchuk & Perlikowski [16] observed that the presence of strongly unstable 
spectrum is not the only possible cause of instability for large N. They observed 
that large numbers of Floquet exponents form bands that have a distance of order 
1/N from the imaginary axis and have a spacing of order 1/N along the imaginary 
axis. In the limit N — > oo these bands form curves after a rescaling of the real 
part by N. The limiting curves, called asymptotic continuous spectrum in [16], were 
defined with the help of a parametric periodic boundary value problem: 

Definition 4 (Asymptotic continuous spectrum). The complex number 7 + icj E C 
(7 € M, u) € [—77, 77)) lies in the asymptotic continuous spectrum, A c , if the 

boundary value problem on [—1, 0] 



has a non-trivial solution y E C 1 ([— 1,0] R n ) for some ip £ R. The quantity ip is 
called the phase corresponding to 7 + iuj. 

Again, the asymptotic continuous spectrum does not depend on N but only on 
A, B and r. 

2.3. A characteristic function for Floquet exponents. We will reduce now 
the study of the various spectra and their relations to each other to a root-finding 
problem of a holomorphic function. The following lemma states the existence of a 
characteristic function for Floquet exponents that at the same time can be used to 
describe the asymptotic spectra A+ and A c . 

Lemma 5 (Characteristic function). There exists a function /i:f!ixfl2CCxCi-> 
C which is holomorphic in both arguments with the following properties: 

1. [i is a Floquet exponent of , i.e. /i E Ejv, */ and only if 



y = [A(t) - fil]y 

1/(0) - 



(ii) 

(12) 



y(t) = [A(t) - iul]y{t) + exp(- 7 - iip)B(t)y((t - r) mod[ _ li0] ) 

y(o) = y(-i) 



(13) 
(14) 



h{n,exp(-{N + r)/i) = 0, 
2. pi is in the instantaneous spectrum T,a if and only if 



(15) 



Mm,o) = 0, 

and, hence, p, is in the strongly unstable asymptotic spectrum A+ if 
h(^,0) = and Re^>0 ; 



(16) 
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3. [i — 7 + ilu is in the asymptotic continuous spectrum A c if, for some phase 

h(iui, exp(— 7 — iip)) = 0. 

The algebraic multiplicity of (i as a Floquet exponent in the statements 1 and 2 
equals its multiplicity as a root in (15) and (16). 

From property 3 the motivation behind the name asymptotic continuous spectrum 
becomes clear: if we have a value 70 + iuj G A c and the corresponding phase </?o, 
and d2h(iuJo, exp(— 70 — i<Po)) is non-zero (which is generically the case) then a whole 
curve 7(0;) +iip(uj) satisfies h(iu>, exp(— 7(0;) — iip(uj))) — for uj sa ujq. These curves 
are the bands of the asymptotic continuous spectrum. Note that from the existence 
of the trivial Floquet exponent [i = we can conclude that /i(0, 1) = 0, which in 
turn implies that 7 = uj = with phase ip — is in A c - In the generic case where 
the trivial exponent is contained in a single curve we call it the critical branch of A c . 

The details of the construction of h, which modifies the general characteristic 
matrices and functions for periodic delay equations from [10, 13], will be given in 
section 3. The general characteristic function constructed by [10, 13] may have poles 
in the complex plane. The modification in Section 3 ensures that these poles of 
h(fi, z) stay in the left half-plane. Hence, the domain VL^ x Jl 2 of h contains all /j, 
and z satisfying Re/i > —R and \z\ < exp(i?) where R > is arbitrary but has to 
be chosen a-priori. Accordingly, statements 1-3 of Lemma 5 are valid only if both 
arguments of h satisfy their respective restriction. However, this is the case in the 
region of interest for stability and bifurcations. 

The introduction of the characteristic function h clarifies how the different spectra 
can be calculated and reduces the analysis of the spectra to a root-finding problem. 
After defining h properly one could even use h to define the corresponding spectra 
by the properties listed in Lemma 5. 

2.4. Main results. With the help of the asymptotic spectra A+ and A c we can 
formulate now a sharp criterion for the exponential orbital stability and instability 
of Mjv, which is our main result: 

Theorem 6 (Stability/Instability). The map (and, hence, the periodic orbit x*) 
is exponentially orbitally stable for all sufficiently large N if all of the following 
conditions hold: 

S-l (No strong instability) all elements of the instantaneous spectrum Y*a have 
negative real part ( this implies in particular that the strongly unstable spectrum 
A+ is empty), and 

S-2 (N on- degeneracy) c^MO, 1) 7^ 0, that is, the Floquet exponent is simple 

for sufficiently large N , and 
S-3 (Weak stability) except for the point \i — with phase tp — the asymptotic 
continuous spectrum A c is contained in {z 6 C : Rez < 0}. 
The map Mn is exponentially unstable for all sufficiently large N if one of 
the following conditions holds 
U-l (Strong instability) the strongly unstable spectrum is non-empty, or 
U-2 (Weak instability) a non-empty subset of the asymptotic continuous spectrum 
A c has positive real part. 

The weak stability condition S-3 is equivalent to stating that for all u) G [— 7r, tt) 
the function z 1— > h^iui, exp(— z)) has no roots with non-negative real part (with 
the exception of z ~ for cj = 0). Similarly, the weak instability condition U-2 is 
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equivalent to stating that h(iu>, exp(— z)) = for some w € [— tt, 7r) and some z € C 
with positive real part. 

Several additional corollaries follow from our analysis: 

Decay rate and dominant frequency. If Mm satisfies S-l-S-3 and is, thus, 
exponentially stable for all sufficiently large N then the decay rate is at most (and 
generically) of order 0(N~ 3 ) and the dominant relative frequency is of order 0(A^ _1 ). 
That is, the dominant non-trivial Floquet exponents are a complex pair of the form 
[-CN- 3 + 0(iV- 4 )] ± i[27r/(JV + t) + 0(N- 2 )}. One branch of the asymptotic 
continuous spectrum touches the imaginary axis in a (generically quadratic) even- 
order tangency at 0, and the dominant non-trivial Floquet exponent lies on this 
branch and next to the tangency. 

Robustness. A map Mm that satisfies S-l-S-3 and is, thus, exponentially stable 
for sufficiently large N according to Theorem 6 remains exponentially stable for all 
sufficiently large N under all perturbations to A, B and r of size less than some e > 0. 
This e does not depend on N because the quantities determining the exponential 
stability do not depend on N but only on A, B and r. This means that periodic 
orbits that are stable for large N are uniformly robust with respect to perturbations 
of the system despite their weak attraction rate of order N~ 3 . 

This permits the conclusion about the coexistence of stable periodic orbits in the 
example in Figure 1. The asymptotic spectra of orbits in the small neighborhood 
IA are small perturbations of the asymptotic spectrum shown in Figure 1(f). In 
particular the curve touching the imaginary axis will also be curved to the left for 
all periodic orbits in U. 

Spectral approximation. Let the instantaneous spectrum have a positive 
distance to the imaginary axis. Then we can guarantee that certain points z in the 
complex plane are in the resolvent set of Mm (that is, the periodic boundary value 
problem (9)-(10) has only the trivial solution for /i = z): 

• every point z in the positive half-plane that is not in the strongly unstable 
spectrum A+ is in the resolvent set of Mm for sufficiently large N. 

• All points of the form j/(N + r) + iuj are in the resolvent set of Mm for 
sufficiently large N if the point z = 7 + ioj is not in the asymptotic continuous 
spectrum A c . 

These two statements about the resolvent set of Mm imply that the spectrum of 
Mm must be close to A+ or (after rescaling) close to A c . The other direction also 
holds if the instantaneous spectrum E^ is not on the imaginary axis: 

• if fj, € A+ has multiplicity k then k Floquet exponents of Mm converge to p, 
for N — > 00 (counting multiplicity). 

• Let 7* + iu* be in the asymptotic continuous spectrum A Cl that is, 

h(iw*, exp(— 7* —itp*)) =0 

for some phase ip* G [— tt,tv). Then we will find Floquet exponents /ijv of Mm 
that satisfy (note that this is small-o) 

RefiN- ^r-f- =o(A^- 1 ), ImMiv-w* = o(l). (17) 
N + T 

Estimate (17) is rather weak. We need and prove a stronger and more detailed 
estimate under the additional condition that c^M^w*, exp(— 7* — 7^ 0. 

Then we have a local root curve 7(w) + i(p(uS) satisfying h(iuj, exp(— 7(0;) — 
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iip(uj))) — for bj near w*, and for sufficiently large N we find algebraically 
simple Floquet exponents fik of Mjv satisfying 



where k are integers such that 2kn/{N + r) is near a;*, and the O-s are smooth 
functions of to. These Floquet exponents y k form a band of discrete complex 
numbers approximating the curve 7(cj) + Hp{uj) of asymptotic continuous 
spectrum. 

Degeneracies. Theorem 6 is sharp except for several degenerate cases that are 
excluded by the conditions S-l-S-3 and U-l-U-2. Degeneracies limiting the region 
of stable periodic orbits are: 

• An element of the instantaneous spectrum has zero real part. 

• The partial derivative c^MO, 1) equals 0. In this case the periodic orbit cannot 
be exponentially stable. 

• (Turing and long wavelength instability) The asymptotically continuous 
spectrum A c touches the imaginary axis in a quadratic tangency at some 
value Hljq. In this case the orbit is not exponentially unstable and it is still 
possible that the periodic orbit is exponentially stable for all TV (this is different 
from the stationary case discussed in [5]) but it may also be weakly stable or 
unstable. A special case is that A c touches the imaginary axis in a quadratic 
tangency in the point y, = with phase <p = n. In this case the periodic orbit 
is still exponentially stable for large N. 

• (Modulational instability) The critical branch of the asymptotic continuous 
spectrum A c , containing the trivial exponent /i = 0, has a tangency with the 
imaginary axis at us = that is of higher order than quadratic. In this case 
the periodic orbit is still exponentially stable according to Theorem 6. 

The characteristic function h constructed in Section 3 is still valid for the degenerate 
cases. However, for a detailed discussion of these degeneracies, one needs not only 
to specify a defining equation for each degeneracy (this is straightforward) but one 
also has to state secondary non- degeneracy conditions. In particular the case of 
instantaneous spectrum with zero real part is somewhat subtle, even though it 
looks similar to the other listed degeneracies of co-dimension 1 at first sight. In the 
analogous situation for the spectrum of equilibria, it turns out that instantaneous 
spectrum iu> with zero real part generically implies a singularity of the asymptotic 
continuous spectrum, 7(^0) = 00, and hence is not part of the stability boundary. 
Due to these reasons, we believe that a comprehensive treatment of the degeneracies 
is beyond the scope of the present paper. 

3. Construction of the characteristic matrix and function. Let us choose a 
constant R > arbitrarily large. We construct a function /ijv(m) that is analytic 
for all fi satisfying Re fi > —R/(N + r), and that satisfies hpf(p) — if and only 
if n is a Floquet exponent of the map Mpf. Hence, this function can be used to 
find all Floquet exponents fi of M/v that satisfy Key > —R/(N + t). Since we are 
interested in the stability of the origin under iterates of Mjy, finding the roots of 
/ijv will then be sufficient. 




(18) 



Re fj, k 



[l + OdN + r)- 1 )} 



7(ImMfc) 

N + T 



(19) 
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We introduce the complex variable ztC and consider the periodic boundary 
value problem for t G [—1,0] 

y(t) = [A(t) - fil]y(t) + z Bit) y{{t - r) mod[ _ li0] ), (20) 
v(P) = y(-l). (21) 

In a first step we will construct a characteristic matrix A(/i, z) for (20)-(21) such that 
the roots of its determinant /i(/i, z) = det A(/n, z) will be precisely those pairs of points 
(/i, z) in some subdomain of C x C for which (20)-(21) has a nontrivial continuously 
differentiable solution y E C 1 ([— 1, 0]; C n ). Thus, by inserting z = exp(— (N + r)fi) 
we will then obtain the characteristic function /iat(/x) in a subdomain of C. 

Consider a partition of the periodicity interval [—1,0] into k intervals of size 1/fc: 



-l + l,-l + 3 -±±) fov./ = () /,■- I. (22; 



Using this partition we formulate a multiple initial value problem (MIVP) for a 
vector of k initial (or restart) values (vq, . . . , Vk-i) T G C nk (similar to multiple 
shooting): 

y(t) = [A(t) - fil}y{t) + zB(t)y({t - r) mod[ _ 1)0] ) (23) 
y (ti)=Vi for i = 0, . . . ,k - 1 (24) 

where t E [—1, 0] and z E C. Notice that in (23) the solution y on the interval [tj, tj + \) 
depends on the solution y(t) in other intervals due to the term y(t — T) mo( j[_ lj o] in 
the right-hand side of (23). 

The main purpose of this partition is to reduce the length of the integration interval 
at the cost of increasing the dimension of the system. Indeed, this construction is 
very similar to the construction in [11], where in the case of rational t a reduction 
to an equivalent system of ODEs could be achieved, and for the case of irrational r 
corresponding rational approximations have been considered. This can be seen by 
introducing Uj(t) := y(t + tj), which satisfy the system of equations 

Uj(t) = [A(t + tj) - [il]uj(t) + zB(t + tj)u mU) ((t + tj - t m{j) - T) mod[ _ lj0] ) (25) 
Uj(0) = Vj (26) 

with j = 0, . . . , k — 1. This system can now be considered as an initial value problem 
on the interval [0, 1/k), and for a solution on this smaller interval [0,1/fc) each 
component Uj(t) represents the solution y(t + tj) in the corresponding subinterval Ij. 
Note that instead of the delayed term Uj((t — T) mod [_ li0 ]) we have inserted a coupling 
to another component u m ^{{t + tj — t m /j\ — i") m od[-i,o]), wnere the index m(j) is 
chosen such that the argument (tj —t m ij\ — T) m od[-i,o] is in the interval [0, 1/k). The 
initial value problem for the coupled system (25)-(26) is an equivalent formulation 
of the multiple initial value problem (23)-(24) and, in a similar way as the boundary 
value problem (20)-(21), contains both delayed and advanced arguments. If r is 
rational then (25)-(26) is a system of ODEs. This was the starting point for the 
rational approximations in [11]. We take another route by showing directly that the 
Picard-Lindelof iteration for (25)-(26) converges, similar to initial- value problems 
of ODEs. In the sequel, we will clarify this point, using the formulation (23)-(24) 
which is more convenient for our purposes. 
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Let us denote by U(t, s, /i) £ R nxn the propagation matrix of the linear ODE 
defining the instantaneous spectrum, (f 1). That is, 

U(t, s, = y(t) where 
y(s) = v and 

y(r) — [A{r) — fil]y(r) for all r £ [s,t]. 
The norm of U (i, s, fi) can be estimated by 

||£/(t,s,M)||oo <exp([||A|U-Re|i](t-*)). (27) 

In (27) we have used the notation \\A\\oo = max te r_ li oi ||oo- We will use the 
same notation for B. In order to clarify in which sense system (23)-(24) is an initial 
value problem and what it means for y to be a solution of system (23)-(24) we 
formulate an integral equation which is equivalent to (23)-(24): 

y(t) = [S(ji)v](t) + z J U(t,s,fi)B(s)y((s-T) modl _ lfi] )ds (28) 

a k (t) 

[S(n)v](t) = U(t,t J , f i)v j Utelj, (29) 
a k (t)=tj ift€lj. (30) 

We note that S(fj,)v and a k are piecewise continuous functions on [—1, 0] ([S (n)v](t) £ 
C" and a k {t) £ R). They are continuous on each sub-interval Ij = [tj, ij+i) but 
have jumps at the times tj. The integral equation (28) is a fixed-point problem for 
y. If we find a fixed point y then y may have discontinuities at the times tj. Thus, 
the appropriate space in which to look for solutions of the fixed point problem (28) 
is the space of piecewise continuous functions with the usual max- norm ||2/||oo : 

Cfc —{y : [—1,0] H> C" : y continuous on each subinterval Ij — [tj,tj + i) (31) 
(J = 0, . . . , k — 1) and lim t y tj y(t) exists for all j = 1 . . . k.} 

The right-hand side of the integral equation (28) is an affine map, mapping C k back 
to itself such that (28) is of the form 

y = S(n)v + zL k (fi)y (32) 

where S(n) : C nk H> C k is defined by (29), and L k {^)y is the linear part of the 
right-hand side in (28) (that is, the integral term). The linear map L k {n) : C k i-> C k 
is continuously differentiable (and, thus, holomorphic) with respect to the complex 
variable fi. A simple estimate for the norm of L k with respect to the || • |joo-norm 
gives us the unique solvability of the fixed point problem (32) (which is actually the 
integral equation (28)-(30)): 

Lemma 7 (Existence and Uniqueness of solutions for IVP). 
Let R > be arbitrary. If we set the number of sub-intervals, k, such that 

k > C{R) :=max{||A|| 00 + R, exp(l + R)} (33) 

then the affine integral equation (28) -(30) has a unique solution y £ C' k for all \x 
and z satisfying 

Refi>-R, |z|<exp(i?) (34) 
and all tuples (uq, . . . ,v k ) T £ C nk . 
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The details of the norm estimate for L k are in Appendix A. The solution y can 
be written as y — [I — zLfc(/i)] _1 S(p)v. When is this solution y (which is only in C'k 
for general v) continuously differentiable on the whole interval? This requirement 
is a linear condition on the tuple v. Let us fix a constant R > and choose the 
number of sub- intervals k > C(R). 

Definition 8 (Characteristic matrix and function). 

We define the characteristic matrix A(fi,z) £ C nfexnfe for the problem (23)-(24) 
as 

v - y(0) 



A(fi,z)v := A(/x,z) 



Vk-l 



vi - lim t /. tl y(t) 



(35) 



y k -i - lim t /. t(c _ 1 y(t) 
and the corresponding characteristic function h(fi, z) as 

h{n, z) := det A(/z, z). (36) 

Note that the integral equation (28) implies that Vj = lirtttvt,. y(t). In this 
sense the values Vj are the initial (or restart) values for the differential equation 
(23). The construction of A gives a well-defined matrix for all z and /x satisfying 
Re/i > —R and \z\ < exp(i?): for a given tuple v we evaluate the unique solution 
y = [I — zL k (fi)]~ 1 S([i)v e Cfc of the integral equation (28) and then we use this 
y to evaluate the right-hand side of (35). The characteristic matrix A is set up 
such that for v in the kernel of A(/x, z) the solution y also satisfies the differential 
equation (23) on the whole interval (and does not have jumps) including the periodic 
boundary conditions: 

Lemma 9 (Differentiability). Let Re /x > — R, \z\ < exp(R) and k > C(R). If the 
tuple v = (vq, . . . , Wfc_i) satisfies 

A(fj,,z)v = (37) 

then y = [I — zL^{^l)]~ 1 S{^,)v, the solution of the integral equation (28), is contin- 
uously differentiable on [—1,0] and satisfies the differential equation (20) with the 
periodic boundary condition (21). Conversely, let y G C 1 ([— 1, 0], C") be a continu- 
ously differentiable solution of (20)-(21). Then the tuple v = (vo, . . . , v^-i) G C fe , 
where v = y(-l), i>i = y(h), v k -i = y{tk-i) satisfies (37). 

Proof. The rows 2 to k of the right-hand side in the definition (35) of A ensure 
continuity of y: since Vj = limi^ tj . y(t) for j = 1, . . . , k — 1 the condition that 
these rows are equal to zero reads limtx^. y(t) = lim f ^ tj y(t), and (by inserting the 
right-hand side of (28)) 

£i'-t-i 



U(t j+1 ,t j ,fi)v i + z / U(t j+1 ,s,fx)B(s)y((s-T) mod[ _ lfi] )ds 



for j = 0, . . . , k — 1. This implies that y is continuous on [—1,0], and that we can 
concatenate all the integral terms in (28) to 



(38) 



y(t) = U(t,-l,n)v + z / U(t,s,n)B(s)y((s-T) mod[ _ lfi] )ds 



The first row of the condition A(/i, z)v = reads Vq = 7/(0) , which makes sure that 
y is periodic at the boundary of [—1,0], and guarantees that the integrand in (38) is 
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continuous at s = r — 1. Consequently, the integrand is continuous everywhere, which 
implies that y is continuously differentiable. Thus, we can differentiate (38) with 
respect to t, which implies that y satisfies the differential equation (20). This, in turn, 
implies that also y{— 1) = y(0) because the right-hand side of (20) is periodic. □ 

Since the map Lk depends analytically on p, the matrix A depends analytically on 
p and z, too (as long as Rep > —R and \z\ < exp(i?)). Calling A the characteristic 
matrix makes sense for the following reason: 

Lemma 10 (Characteristic function for (6) [10]). Let N > 1 and R > be given, 
and choose the number of sub-intervals, k, greater than C(R). Then for all p 
satisfying 



the following equivalence holds: p is a Floquet exponent of the time-1 map Mm of 
the periodic linear DDE (6) if and only if 



h N (p) := hiii, exp(-(N + r)fi)) = det A(p, exp(-(iV + r)p)) = 0. (39) 



The algebraic multiplicity of p as a Floquet exponent of Mm is equal to the order of 
as a root of hpf. 

Proof. That p is a Floquet exponent of (6) if and only if A(^,exp(— (N + t)/x)) 
has a non-trivial kernel is clear because Floquet exponents of Mm were defined 
as those complex numbers for which (20)-(21) has a non-trivial solution. The 
statement about the multiplicity of pi follows from arguments similar to [3], which 
are laid out in detail in [10] (an equivalence between the eigenvalue problem for 
Mm, exp(/i)u — Mmv — 0, and the algebraic equation A(p, exp( — (N + r)fi))v = 
is constructed in [10]). □ 

Both functions, h and hjy, are real analytic (that is, their expansion coefficients 
are real numbers since A, B and r are real). 

For any given R > we have chosen k > and constructed a characteristic 
function h^ifi), defined on the half plane 



the roots of which arc precisely the Floquet exponents (counting multiplicities) of 
the period map Mm for the linearized delay differential equation (6). 

As we are interested in the location of Floquet exponents of the time-1 map Mm 
to the right or close to the imaginary axis we have reduced the eigenvalue problem 
to a study of the asymptotic behavior of roots of the holomorphic function hM for 
N — > oo. Furthermore, the asymptotic spectra are defined as roots and root curves 
of h by construction of ft.: p > is in the strongly unstable spectrum A.+ if and only 
if h(p, 0) = 0, and p = 7 + iuj is in the asymptotic continuous spectrum A c if and 
only if h(iui, exp(— 7 — iip)) = for some phase tp € R. 

This proves Lemma 5 and allows us to follow an approach similar to [17, 5] in 
the following sections. Note that in the case where the period and the delay are 
rationally dependent, the existence of a suitable characteristic function follows much 
easier (see [11]), but the structure of the asymptotic spectra discussed below will be 
not affected. 
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4. Strongly unstable spectrum. One immediate consequence of Lemma 5 is 
the statement asserted in [16] for the relation between eigenvalues of Mjv and the 
strongly unstable spectrum A+. Note that /z is in the strongly unstable spectrum 
A+ if and only if h([i, 0) = 0: if one sets z = in (20) the differential equation 
reduces to (11), the expression defining A+. For large N unstable Floquet exponents 
of Af/v either approach the imaginary axis or A+: 

Lemma 11 (Convergence to strongly unstable spectrum). 

Let Re > 0. If /i(/i0i 0) 7^ then there exists a No such that /j,q lies in the resolvent 
set of the time-1 map Mn for all N > Nq . 

If /io is a root of multiplicity k of the function h{-, 0) then every sufficiently small 
neighborhood U of fXo contains exactly k Floquet exponents (counting multiplicity) of 
Mjy for all N > N (N depends on U). 

The resolvent set is the set of all complex numbers fi for which the map exp(/i)/ — 
Mm is an isomorphism (characterized by /ijv(m) 7^ for Re/i > —R/(N + r)). 

Proof. The statement for h(fio, 0)^0 follows from expansion of h(fj,o, •) in z = 
since z — exp(— (iV + r)fi ) — > for N — > 00 if Re [i > 0. 

Let fj,o be a root of h(-, 0) of multiplicity k. For any sufficiently small S > 0, /io is 
the only root of h(-,0) inside the ball Bs(no) of radius S around fiQ. In particular, 
/&(//, 0) ^ on the boundary of Bs(fio)- Then h(n, exp(—(N + r)fi)) is also non-zero 
on the boundary of Bs(no) for sufficiently large AT, and the logarithmic derivative of 
h{(i, exp(— (N + t)h)) converges for N —> 00: 

^[fe( M ,exp(-(jV + r)^))] ^ 

h{ii,exp{-{N + T)n)) 
dih(n, exp(-(iV + t)^)) -(N + t) cxp(-(7V + T)n)d 2 h(ii, expj-{N + r)/z)) 

ft(/i,exp(-(JV + T)M)) 

^tv^oo ^tt^-t^ if Re/i > for all /u G B s (n ). 
h{n, 0) 

Thus, the line integral of the logarithmic derivative of h(/i, exp(— (JV+r)/x)) along the 
boundary of Bs(hq), which counts the roots (and poles), also converges for N — > 00 
to the logarithmic derivative of h(-,0). Since the line integral of the logarithmic 
derivative is an integer it must be constant for N — > 00, and, hence, /ijv(/^) must 
have the same number of roots (counting multiplicity) inside Bs(fio) a s h{-, 0), which 
is k. □ 

Lemma 11 shows that is exponentially unstable for all sufficiently large N if 
the strongly unstable spectrum is non-empty, which is condition U-l of Theorem 6. 

5. Asymptotic continuous spectrum. Due to Lemma 5 the asymptotic contin- 
uous spectrum is given as those fi = 7 + iu> for which one can find a phase (p G K 
such that h(iw, exp(— 7 — iip)) — 0. 

Thus, we expect the asymptotic continuous spectrum to come in curves: Let 
70 +iuio G A c be a point in the asymptotic continuous spectrum and let ipo G [—it, tt) 
be its phase. By definition ^ — 7o + ifo € C satisfies = h{iu)Q, exp(— ^q))- If 
d2h{iu>o, exp(— /Xo)) ^ then there is a root curve (i(u>) of complex numbers satisfying 
h(iu, exp(— fJ-(u)))) = going through jito- Hence locally (for u in a neighborhood of 
wq) there is a curve u) > (Re (j,(u>) + ioS) G A c through 70 + iioo for to near wq • 
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We assume for the remainder of the section that the instantaneous spectrum 
has a positive distance to the imaginary axis. The idea behind the construction of 
the asymptotic continuous spectrum is that for Floquet exponents fi close to the 
imaginary axis, that is, for /i of the form 

*=jfc+iu> (40) 

with a bounded factor 7 in the real part, the roots of the characteristic function 
h>N{fi>) = h(/i, exp(— (N + t)h)) converge to a regular limit for N — > 00 after inserting 
the scaling (40) . This can be made more specific: if the instantaneous spectrum 
Tia has a positive distance from the imaginary axis then all Floquet exponents that 
are not converging to the strongly unstable spectrum have a real part less than 
Rz/iN + t) where R2 does not depend on N: 

Lemma 12 (Convergence to the imaginary axis). Assume that the elements of the 
instantaneous spectrum have a positive distance 2e > from the imaginary axis. 
Let us denote the points in the strongly unstable spectrum, A.+ , by . . , pLj. There 
exists a constant R2 > such that all Floquet exponents fi of Mjv satisfy either 
/1 € B e (jJ,j) for some j, or 

Re ^<prtr (41) 

for allN > 1. 

Proof. For all complex numbers fi that have non-negative real part and are outside of 
the balls B e (p,j) the matrix I — U(0, — 1, /x) is invertible and the inverse has a uniform 
upper bound (remember that U is the propagation matrix of y — [A(t) — jJLl]y): 

IKi-f/Co,-!,//)]- 1 !! <c. 

We also know that n is a Floquet exponent if the integral equation (38) (which is 
equivalent to the differential equation (37)) has a non-trivial periodic solution y(t) 
for z — exp(— (N + r)/x). This implies (y(— 1) = y(0) = v n ) 

v = U(0, -l,fi)v + z J U(0, s,n)B(s)y((s - r) mod[ _i, ])ds, 
and, hence, 

r ,0 
y(t) =z [I - U(0,-l, nT 1 / U(0, s, n)B(s)y((8 - r) raod[ _ 1> o ] )ds+ 

(42) 

U(t,s,ix)B(s)y((s - r) mod[ _ l!0] )ds 

1 -1 

The factor of z in the right-hand side of this fixed-point problem for y is a linear 
operator K (fi) on the space of continuous functions. K{fx) is uniformly bounded 
for all pi that have non- negative real part and are outside of the balls B e (fj,j). Let 
us denote an upper bound on the norm of K(fj.) by C\. Then for all z satisfying 
\z\ < C-f the integral equation (38) cannot have a non-trivial solution. Since, for 
Floquet exponents of Mjsr, z equals exp(— (N + t)/i) this means that fi has to satisfy 

|exp(-(iV + r)/i)| >Cr\ 

and, thus 

Reu < — 

r ~ N + T 
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for \z\ to be larger than C-f 1 . Consequently, if we choose R 2 > logCj then a /i that 
has non- negative real part and lies outside of the balls B e (fij) must satisfy (41) in 
order to be a Floquet exponent of Mm- □ 

A corollary of the construction described by equation (42) is that for each ui 
the intersection of the asymptotic continuous spectrum A c with the horizontal 
line {/i G C : Im^ = lo} can be described as a set of eigenvalues of a compact 
operator K(iuj). Thus, if ReXU ^ the asymptotic continuous spectrum forms 
curves parametrizable by w, which may only have poles at — oo: 

Corollary 13 (Asymptotic continuous spectrum consists of curves). Assume that 
the instantaneous spectrum Y>a has a positive distance to the imaginary axis. A 
point fx is in the asymptotic continuous spectrum if exp(fi) is an eigenvalue of the 
compact linear operator K(iuj) : C([— l,0];C n ) <-> C([— 1,0]; C") given by 



K(\u)y :-- 



[I -[/((), -l,^)]" 1 / U(0,s,icj)B(s)y((s-T) mod[ _ 1<0] )ds+ 



U(t, s, iu)B(s)y((s - r) mod [_i )0 ])ds 



In particular there exists a constant R3 > so that the asymptotic continuous 
spectrum lies to the left of the vertical line {zeC: Rez = -R3}. 

The matrix U in the definition of K is the monodromy matrix of the instantaneous 
problem x(t) = [A{t) — iuj]x(t). The inverse of / — U(Q, — l,isJ) exists because of 
our assumption that the instantaneous spectrum, S^, does not contain points on 
the imaginary axis. The corollary implies that the asymptotic continuous spectrum 
consists of curves even in points where the regularity condition on d 2 h 7^ is violated. 
In these points K(iuj) has eigenvalues of finite multiplicity larger than 1 such that 
finitely many curves of the asymptotic continuous spectrum cross each other (or are 
on top of each other). 

All Floquet exponents of Mm that do not converge to the strongly unstable 
spectrum for N — > 00 (these can be at most n, equaling the dimension of A(t)) are 
either stable, or they satisfy restriction (41) on the upper bound. Thus, apart from 
the strongly unstable spectrum, the multipliers, which are of interest from the point 
of view of stability and bifurcations, lie in the strip 

C N ■■= \n € C : -— — < Rc^i < h - 



N + t ~ N 



and have the form 



7 

/it = — h ioJ where 7 6 \—R, R2] and w G — 7r, tt) (43) 

Jy + r 

(after shifting them into the strip {z : Imz € [— 7r,7r)} by subtraction of an integer 
multiple of 27ri). As we focus our discussion on Floquet exponents of the scale (43) 
from now on it makes sense to introduce the notation of a scaled Floquet exponent. 

Definition 14 (Scaled Floquet exponent & resolvent set). 

A complex number fi = 7 + \uj with 7 6 [— R, R 2 ] and u € [— 7r, ir) is called a scaled 
Floquet exponent of if j/(N+t) +iu> is a Floquet exponent of Mjy. Similarly, 
H is in the scaled resolvent set of Mm if 7/ (N + t) + iu is in the resolvent set of 
M N . 
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We can now formulate precisely in which sense the asymptotic continuous spectrum 
A c is the limit of the spectra of Mjy. First, we make a statement about resolvent 
sets. 

Lemma 15 (Points distant from A c ). Let the instantaneous spectrum Ha have 
a positive distance to the imaginary axis. If 70 > —R and 70 + iujQ is not in the 
asymptotic continuous spectrum A c then 70 + iujQ is in the scaled resolvent set of 
Mpj for sufficiently large N . 

Proof. Lemma 15 is a simple consequence of the fact that, if h(iu}Q, exp(—jQ—iip)) 
for all phases tp G [— tt, n] (which defines points 70 + iujQ (£ A c if 70 > —R) then 

h N ~ + iuoj = h (ioJo + ^ , exp (-70 - i(N + t)uj )^ ^ 

for all sufficiently large N. Lemma 5 implies that 70/ (-/V + t) + iojQ is in the resolvent 
set of Mm for all sufficiently large N. □ 

Lemma 15 implies that all spectrum of Mn near the imaginary axis has to be 
close to A c even after blow-up of the real part. The other direction, that every point 
of A c in the strip —R < Re /1 < R2 is approached by scaled Floquet exponents of 
Mjv, is also true (again under the assumption that the instantaneous spectrum £4 
is not on the imaginary axis). We prove this direction first for regular root curves 
of h because this is the most common case, and the missing piece for the stability 
criterion. Our estimate is slightly sharper than mere approximation to make it useful 
for our proof of the stability criterion. 

In short, Lemma 16 below states that any regular curve of asymptotic continuous 
spectrum, 7(0;) + iui (with its corresponding phase tp(u>)), is approximated by scaled 
Floquet exponents of the form 7^. + iw^ where 

7fc = [I + O^ + t)- 1 )]^), (45) 

k are integers such that 2kn/ (N + r) is near to, and the O-s are smooth functions of 
lj. This proves the stronger estimate (18)— (19), given in the non-technical overview 
in Section 2. 

Lemma 16 (Convergence to regular curves of A c )- 
Assume that the triplet (w* , 7* , if* ) satisfies 

h(iu>*, exp(— 7* — iip*)) = 0, 92/1(1^*, exp(— 7* — iip*)) ^ 0. 

Let be sufficiently large and S > be sufficiently small. We denote the unique 
regular root curve of ^ 1— > h(iw, exp(— fx)) for u near through (u) = = 
7* + itp*) by 

£l(u>) = 7(0;) + itp(<jj). 

For all N > there are scaled Floquet exponents /J,jsr,k of Mn near 7* + icu* . All 
scaled Floquet exponents /ijv.fe that are sufficiently close to 7* + iw* are algebraically 
simple and have the form 

f-N.k = lN{^k) + iuk, (46) 

where: 
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k is any integer satisfying 
2kn 



e (w» - 6,u* + S), 



N + T 

• cjfc is £/ie unique solution of the fixed point problem for lo 

^) (47) 
the functions 7/v and <f>N are perturbations 0/7(0;) and (p{oS) of the form 



7Ar(oj) 



1 / 1 



7(w) (48) 



<£wM = <£(<*>) + ^ - Imp ^cj, 7(0;) (49) 

where g(u,e) is a smooth complex-valued function, which is independent of N 
and defined for lo G (w* — S, lo* + (5) erne? e G [0, 1/N*). 

We notice that the scaled Floquet exponents of Afjv he on bands: they are on the 
curve given by 7/v(u;) 4- ilo, which is a small perturbation of the curve of asymptotic 
continuous spectrum 7(0;) + ilo. The spacing between scaled Floquet exponents along 
this band is given by the fixed point equation (47). The fixed point problem (47) is 
only weakly implicit since the right-hand side terms containing u) all have a pre-factor 
1/(N + r), which is small. Hence, expression (45) for 7^ follows immediately from 

(48) and expression (44) for lo^ follows from (47) and (49). 

Proof, (of Lemma 16) We know that j(u>) + i</?(w) is a regular root curve of 

h(iw, exp(— 7 — 'up)) = 0. (50) 

The root problem 

h (iu) + 67, exp (—7 — 'up)) = 

is for small e a small (order e) perturbation of the root problem (50). Thus, for 
sufficiently small e, a root curve of the form %(u)) + \(p t (uj) exists for lo in some 
neighborhood of lo* , and it has the form 

%(u) + itp e (co) = [1 + eg (lo, e)] j(lo) + i<p(w), (51) 

where g(u>,e) is a smooth complex- valued function defined for w in a small neighbor- 
hood of and e € [0, e max ) (with some e max > 0). Note that the error term contains 
a factor j(lo), making the error equal to zero on the imaginary axis. (See Appendix B 
for details of how to extract this factor from the error.) Inserting e = (N + t) _1 
and labeling the curves 7 e (w) as 7jv(w) and tp e as (p^ gives the definitions (48) and 

(49) in the lemma. Correspondingly, we make an initial choice for the minimal AT, 
, as 1 /e ma x ■ 

A point on the curve 7at(cj) + ilu is a scaled Floquet exponent of Mm if and only 
if its imaginary part 10 satisfies 

exp(— \lo(N + t)) = cxp(— itfN ((*;)) and, thus, 

(N + t)lo — ipn(lo) + 2kir for some integer k E Z. 

After dividing by A^ + r this becomes the fixed point equation (47) for lo. For which 
k does this fixed point problem have a unique solution? 
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We choose a neighborhood U of w* of the form (ui* — 26, + 26) such that cp^ 
is well-defined for all u> G U and all N > N*, and satisfies < 1/ for some 

constant L (which is independent of N). Next, we increase iV* such that 

L < 1 for all N > iV*, and 

< 6 for all N > TV* and w G W. 



7V + T 



iV- 

Then the right-hand side of the fixed point problem (47) is contracting with a 
Lipschitz constant L/(N + r) for all w G W and it is mapping U = (w* — 2(5, a;* + 2<5) 
back into itself: 



<p N ((j) 



2kw 



N + T N + T 

for all u) G U and all AT > AL if 



< 



2fc?r 



<pjv(w) 



N - 



2kir 



N 



< 26 



N - 



< 6. 



(52) 



Thus, the Banach Contraction Mapping Principle guarantees that (47) has a unique 
solution ujk for all integers k satisfying (52). 

Finally, we confirm the algebraic simplicity of the scaled Floquet exponents 
HN,k — TivC^fc) + luJ k by checking the multiplicity of the corresponding root of h^: 
the derivative of h^, divided by N + r, is 

h! N (z) d\h(z, exp(— (N + t)z)) , ... . , „ , . . . ... 

= y + V r — - exp(-(7V + r)z)d 2 h(z, exp(-(7V + r)z)). 

Inserting 7jv(wfc)/(A r +r)+ia;fc for 2 on the right, and the relation exp(— i(N+r)uj k ) = 
exp(-iip N (uJk)) we get 

' -h' N (z)=0(N- 1 )- 



N - 

-exp(-7Ar(w fe ) - iip N (uj k ))d2h [ iuo k + ^ k , exp(— Tjv(^fc) ~ i<fiN(^k)) 

\ iV + T 

=0(A r_1 ) - exp(-7Ar(a;fc) - i(p N {Lo k ))d 2 h(iuj k , exp(-7Ar(w fc ) - i(p N (uJ k )))- 

Since is in the neighborhood W of w* in which is non-zero the overall derivative 
is non-zero for sufficiently large N . □ 

Lemma 16 is based on a perturbation argument assuming that the complex 
function h(iw, exp(— 7 — i<p)) has a regular root curve cj h-> j(u>) + \(p{oS). Thus, it is 
also valid if the instantaneous spectrum, S^, does not have a positive distance to 
the imaginary axis as long as one restricts consideration to Floquet exponents of the 
form 7/ (N + t) + iu> (with bounded 7). Positive distance of to the imaginary 
axis merely ensures that all Floquet exponents [i with Re /1 > —R/(N + r) are of 
this form (except for those approximating the strongly unstable spectrum A+). 

Lemma 15 and Lemma 16 about the approximation of the asymptotic continuous 
spectrum, together with Lemma 11 about the approximation of the strongly unstable 
spectrum, are the tools that we need to prove the criterion for asymptotic stability 
from Theorem 6. Before turning to asymptotic stability let us prove the remaining 
statement about spectral approximation of the asymptotic continuous spectrum Ac- 
Let 7* + icj* be an element of A c with phase 93* (that is, h(iu)*, exp(— 7* — i<p*)) = 0). 
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Then we find scaled Floquet exponents of Mjq that approximate 7* + iu;* even if 
the non-degeneracy condition <?2/i(iw*, exp(— 7* — iy*)) ^ is not satisfied: 

Lemma 17 (Approximation of A c ). Assume that the instantaneous spectrum Y,a is 
not on the imaginary axis. Let 7* G [— R, R2] and let 7* + iw* be in the asymptotic 
continuous spectrum A c . Then there exists a sequence "fpj + iuipj of scaled Floquet 
exponents of such that 

7jv + i^AT — > 7* + iw*. as N — > 00. 

Lemma 17 covers the claim of Theorem 6 about exponential instability of M/v 
for sufficiently large N under the condition of weak instability U-2. 

Proof. Since the instantaneous spectrum Ha does not contain points on the imaginary 
axis, exp(7 >t + iip*) is a non-zero eigenvalue of the compact operator K{iu*) as 
introduced in Corollary 13. Consequently, exp(7* + i^*) is isolated and has finite 
multiplicity, which implies that 7* + iip* has finite multiplicity as a root of z ^ 
h(iu*, cxp(— z)). Let us define k(N) for large N as 

2kir 

k(N) = greatest integer k such that < w*. 

v / f N + t - 

By construction of k(N) we have that 

27rfc(iV) 
hm = w* . 

N->oo N + T 

Since z h-> h(iw*, exp(— z)) has an isolated root at 7* + i</?*, the functions 

, f2irik(N) z . . 

z ^H^T^ + ^T7' cxp( ~ z) 

also have roots z^ = "fN + i^TV for sufficiently large N which converge to 7* + itp* 
for — > 00. Let us define 

2nk(N) + ip N 

N + T 

Then, by construction, 7^ + iwjv is a scaled Floquet exponent of Mjv since 
' "f N + 2-Kik{N) + i(p N 



v _ _ -,exp(-7jv - (2nk{N) + tp N )i) J = 0. 

Moreover, 77V — > 7* and ujn — > oj*, which proves the claim of the lemma. □ 



6. Asymptotic stability for large delay. The convergence results for the strongly 
unstable spectrum in Lemma 11 and the asymptotic continuous spectrum in 
Lemma 15 and Lemma 16 can be combined to give a criterion for the stability 
of the periodic orbit of the original nonlinear system (5) depending on the triplet 
(A, B, t) that ensures stability of x* for all sufficiently large N . If the instantaneous 
spectrum T,a has a positive distance to the imaginary axis then the point is part 
of the asymptotic continuous spectrum because is a Floquet exponent for all N . If 
92/i(0, 1)^0 then a regular curve 7(0;) +icu of the asymptotic continuous spectrum is 
passing through (that is, 7(0) = 0). This curve at least touches the imaginary axis 
because 7(0) = and 7'(0) = is an even function, thus, all odd derivatives 

of 7 are zero). 
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Lemma 18 (Asymptotic stability). Let the triplet (A, B, r) be such that all elements 
of its instantaneous spectrum have negative real part. Furthermore, we assume 
that d2h(0, 1) ^ 0, and that for the asymptotic continuous spectrum (including the 
corresponding phase ip ) 

A ClV := {((J,j,(p) : h(iu>, exp(~7 - iip)) = 0, 

u G [-7T,7r), 7 G [~R,R 2 ], ip G [-7T,7r]} 

the point (7,0;, (/j) = (0,0,0) is the only point with 7 > 0. Then the map Mjy is 
orbitally exponentially stable for all sufficiently large N . 

Note that the assumptions of Lemma 18 exclude the case h(0, —1) = since this 
would mean that 7 = u> — 0, tp — n is in A Ctip . 

Proof. (Lemma 18) Since c^MO, 1) ^ we know that the Floquet exponent is 
simple for Mn if N is sufficiently large. Also, since the instantaneous spectrum is 
in the negative half-plane the strongly unstable spectrum is empty. Hence, Mjv is 
the return map of a stable periodic orbit x* if it has no non-zero scaled Floquet 
exponent 7 + ius for which 7 G [0, R 2 ] (and ui G [— n, ir)). 

Proving the statement by contradiction, we assume that, for a sequence of 
increasing AT, Mn has a scaled Floquet exponent 7jv + iwjv ^ where jn G [0, R?]. 

The sequences ( , Jni^n, ( Pn) where (pN = (N + t)ujn mod [— 7r,7r) must have 
accumulation points. Without loss of generality we pick our sequence such that it 
converges to one of these accumulation points, say (7*, o>», (p*). Since hjv(7jy/(iV + 
r) + iujn) = we have by definition of ft, at and <pn'- 

= hi icj _/y + ^ — , exp(— 7jv — i<Pn) ) , and, thus, by continuity of h 
\ N + t J 

= h(iuj*, exp(— 7* — i<p*)). 

Consequently, the accumulation point must be an element of the asymptotic con- 
tinuous spectrum A c , v . Since 7 at > for all N of the sequence, 7* must be greater 
or equal 0, too. By assumption, the only element of A c>v with non-negative 7 is 
7* + ioj* = 0. Thus, 7* = = = 0. 

We choose e > sufficiently small such that we can apply Lemma 16 to (7* , tp*, w*) . 
We have that ojn G (— e, e) and 7^ + iip N G B e (0) for sufficiently large N of the 
sequence (-f N ,uj N ,tp N ) (how large N has to be, depends on e). Since d 2 h(0, 1) 7^ 
this guarantees that 7 at lies on the curve 7 at (w) given in (18) in Lemma 16: 



1 / 1 



JV + t \ N 



Since for all sufficiently large N the factor in front of 7(0; jy) is positive, jn^n) must 
have the same sign as j(u/n), which is negative if ujn ^ due to the assumptions of 
the lemma. Since 77V is assumed to be non-negative, this implies that wjv = "fN = 0. 
Thus, the scaled Floquet exponents 7jv + i^jv of Mpj are zero for the converging 
sub-sequence, which is in contradiction to our assumption 7^ + iu>N =/= 0. □ 

Lemma 18 proves the exponential stability claim of Theorem 6. 

7. Conclusions. We have shown that Floquet exponents of periodic solutions of 
delay differential equations (1) with large delay can be approximated by a set of 
continuous curves (asymptotic continuous spectrum) that are independent of the 
delay, and a finite set Floquet exponents (strongly unstable spectrum) . Although 
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the structure of the spectrum is shown to be similar to the case of equilibria [5] , 
there are some unique features, which occur specifically for periodic orbits. Our 
results are based on the construction of the characteristic function, the roots of 
which give Floquet multipliers of the periodic orbit. 

Using the asymptotic spectra we have been able to provide necessary and sufficient 
conditions for the exponential stability of periodic solutions for all sufficiently large 
delays. Our results are applicable to the case when the delay r is large compared to 
the period T of the solution. In this case, the large parameter N, which controls 
precision of the asymptotic approximation is proportional to r/T. 

Let us mention some of the specific features of the spectrum. In contrast to 
the equilibrium case, the asymptotic continuous spectrum of Floquet exponents for 
periodic solutions contains generically a curve with a tangency to the imaginary 
axis (see Figure l(e,f)). We have proved that even in the presence of this tangency, 
the stability (or instability) of the asymptotic continuous spectrum implies the 
exponential stability (resp. instability) of the corresponding periodic orbit. We 
have shown that the generic decay rate of perturbations of the exponentially stable 
periodic orbit of system (1) is of the order JV~ 3 . 

From the practical point of view, our results can be useful for studying periodic 
regimes in applications that involve feedback with large delays, for example, semi- 
conductor lasers with optical feedback [4, 17], or systems with feedback control 
[9]. 

From a mathematical point of view our result may provide a rigorous approach 
to proving the existence of a large number of stable rapidly oscillating periodic 
solutions for some special cases in which the asymptotic spectra can be computed 
explicitly. This would provide a contrast to the results for scalar feedback equations 
[14]. 
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Appendix A. Proof of Lemma 7. We have to estimate the norm of L k with 
respect to the || • norm. The operator L k was defined as 

[L k {p)y] (t) = \ U(t,s,fx)B(s)y((s-T) mod [- lfi ])ds 

Ja k (t) 

mapping a piecewise continuous function y G C k back into C k . Using the norm 
estimate (27) we can estimate the norm of L k by 



\Lk(u)v\\co < max 



f ||[/(m,m)IU|£|UMI 

Ja k (t) 



< max / eM[\moo~'Refj](t-s))ds-\\B\\ 00 \\y\\ 00 . (53) 
te[-i,o] Ja k (t) 

We distinguish two sub-cases depending on the sign of ||^4||oo — Re/x: 

Case 1. If || Alloc < Re ft then the integrand in (53) is bounded by unity such that 

||ifc(A*)||oo < (* - afc(i))||5||co < £||B||oo, (54) 

taking into account that the length of the integration interval [afe(i),t] in the right- 
hand side of (28) is less than 1/k for all t e [—1,0] by construction of ak(t) (see 
(30)). 

Case 2. For the case Reft, < ||^4||oo we can bound the whole integral by (note that 
for any a > 0, t > the inequality [exp(ai) — l]/o < atexp(at) holds) 

cxp([||A|| 00 -Rcfi] (t-a k (t)))-l < 
Halloo -Re/x 

(t - a k (t)) exp ([HAIloo - Re^] (t - a k (t)) . (55) 

One of the conditions of Lemma 7 was that k > ||A-||oo + R- Thus, if Re fi > — R we 
have that 

[WAW^-Re^it-a^t)) < [WM^ + R] (t - a k {t j) < - [WA^+R] < 1 (56) 

since < t — a k (t) < 1/k by definition of a k . Inserting (56) into (55) the integral 
term in (53) is bounded by exp(l)/fc such that 

\\L k (n)\\oo < ^^HSlloo (57) 
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Inserting (57), the more pessimistic of the two estimates (54) and (57) for both 
cases, into the upper bound for Lk we get 

||£ fc (A»)||oo<exp(l)H^. 

Condition (33) on k from Lemma 7 (requiring that k > ||B||oo exp(l + R)) implies 
that the norm of zL^ is less than 1. Consequently, I — zLk(fJ-) is invertible such that 
the fixed-point problem (32) has a unique solution for all tuples v. □ 

Appendix B. Multiplicative perturbations. In Lemma 16 we had a triplet 

(cj ,7o, (fo) such that 

h(iu) , exp(-7o - iip )) = 0, and d 2 h(iu , exp(~7 - iip )) ^ 

(h was an analytic complex function in both arguments). One has a regular local 
curve "/(uj) + \tp(ui) of complex numbers near 70 + iujq satisfying 

h(iw, exp(— j(oj) — iip(uj))) = 
for all uj ps wo- Then, Lemma 16 claims, the regular root curve 7 e (w) + i<p e (uj) of 

h(h) + ej, exp(-7 - itp)) = 0, (58) 
which exists for small e, has the form 

j e (u) + \<p e {u) = 7M + i ( y 5 ( w ) + e 3e( w )7( w )- (59) 

The emphasis in (59) is on the factor 7(0;) in the error term eg e (u>)~/(uj), which comes 
from the special type of perturbation in equation (58) defining the curve. Note that 
g e {u>) is complex. 

This fact is a special case of the following general statement: 

Lemma 19. Let e be small, the function f : E™ x M. m i-> K m be smooth, /(a;*, J/*) = 
and $2/(2:*, V*) be invertible. Let A G R nxm be a matrix. Then the curve y t {x) 
defined implicitly by 

f(x + eAy e ,y e ) = (60) 

for x w has the form 

y t (x) = [I + eg(x,e)A}y (x) (61) 
where g is a m x n matrix depending smoothly on x and e. 

Note that yo(x) is the curve defined implicitly by f{x,yo) = (putting e = in 
(60)). 

Proof. The Implicit Function Theorem guarantees that y e {x) exists for small e and 
x fa x*, and that it has the form y e (x) = yo(x)+eh(x, e). Subtracting the expressions 
f(x + eAy e , y e ) and f(x, yo(x)), which are both zero, from each other, and applying 
the mean value theorem we obtain 

= eD x {x, e) Ay e + D 2 {x, e)[y t - y (x)] (62) 

where D\ and D2 are the averaged derivatives: 



Di = I d 1 f(x + seA(y (x) + eh(x,e)),y (x) + seh(x,e))ds 
Jo 

D 2 = I d 2 f(x + seA(y (x) + eh(x,e)),y (x) + seh(x,e))ds 
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Note that we have replaced y e by yo + eh(x, e) inside the arguments of D\ and D 2 - 
Since D 2 is invertible for small e and x ~ re* we can rearrange (62) for y e (dropping 
the arguments x and e from D\ and D 2 ' 



Note that (63) is not an explicit definition of y e (x) but rather a fixed point 
problem for y e because y e occurs on the right-hand side as well (via the unknown 
function h). However, the Banach contraction Mapping Principle can be applied to 
the fixed point problem (63) to produce an explicit definition of y e . 

If we treat the complex numbers as a two-dimensional vector space then oj (treated 
as complex number) plays the role of x € K 2 , 7 + i(p € C = M 2 plays the role of y, 
(x, y) h-> h(\x, exp(— y)) plays the role of /, and y i-> — i Re y plays the role of A. 
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E-mail address: matthias.wolfrumOwias-berliii.de 
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y e = [1 + eD^DxA] y (x) 




(63) 



which is of the form (61) as claimed by the lemma. 



□ 



